Numerical Solution of Hard-Core Mixtures 
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We study the equilibrium phase diagram of binary mixtures of hard spheres as well as of parallel 
hard cubes. A superior cluster algorithm allows us to establish and to access the demixed phase for 
both systems and to investigate the subtle interplay between short-range depletion and long-range 
demixing. 



PACS numbers: 64.75. +g 61 .20.Gy 

Liquid binary mixtures pose an important, yet easily 
formulated problem: Imagine a 3-dimensional box of size 
L (volume V = L 3 ) occupied by objects of two differ- 
ent types a, b with packing fractions r\ a and rjb- Will the 
system remain homogeneously mixed or will it phase- 
separate? 

This demixing problem can be readily analyzed when- 
ever there is an obvious free energy imbalance between 
the homogeneous and the phase-separated system. In 
many cases of practical importance, such an imbalance 
is due to the interactions, for example caused by elec- 
trostatic screening. Even in the absence of interactions 
(other than by a hard-core term), a free energy differ- 
ence can be caused by an entropic contribution. In "non- 
additive" mixtures, the, say, demixed phase may be able 
to pack space more densely. At high overall packing frac- 
tions, the system will then be phase-separated. 

Systems of impenetrable large and small spheres or 
cubes belong to the class of additive mixtures. In these 
systems, the distance r a b of closest approach between, 
say, two spheres of radii r a and r& satisfies: r Q & = r a + rb, 
so that the abovementioned simple entropic effect is ab- 
sent. Nevertheless, it has been understood for a long 
time that even additive mixtures are subject to an en- 
tropic "depletion force" Q: two large particles may ap- 
proach sufficiently closely for the small ones to be ex- 
pelled from the interspace between them. In that sit- 
uation, an osmotic (partial) pressure difference (of the 
small particles) builds up and pulls the big particles even 
closer together. The depletion force is strongly attractive 
at very short distances r between the large particles (for 
2 r„ < r < 2 r a + r&), but finally turns out to be quite 
long-range in nature 0] ||. Approaches to integrate out 
the small particles are thus problematic. As a prototype 
of additive mixtures, it is thus of great interest to com- 
pletely analyze the microscopic model of large and small 
spheres or cubes and to understand whether it will even- 
tually lead to phase separation. Precisely this question 
has remained hotly debated even in recent years. 

The theoretical framework for the mixture problem has 
traditionally relied on the solutions of 'closure approx- 
imations' to the Ornstein-Zernicke integral equations. 



Classic work of Lebowitz and Rowlinson Q, performed 
more than 30 years ago, first showed that the Percus- 
Yevick closure of the integral equations did not lead to 
phase separation. This exact statement was thought to 
reflect the true behavior of the system until Biben and 
Hansen jj| challenged the view by showing that an insta- 
bility of the mixture was predicted by a different choice of 
the closure. Since the closure approximations are contra- 
dictory (and fundamentally uncontrolled), it is of prime 
importance to resort to independent checks. However, 
numerical simulations have been notoriously difficult, es- 
pecially for objects very different in size. Numerical evi- 
dence for a phase transition has, to our knowledge, only 
been obtained in a lattice model of hard cubes || . 

The situation has thus been very confusing. One is lead 
to agree with the author of ref. 0: the field would ben- 
efit from an "Ising model of liquids" , an exactly solvable 
full-complexity model against which the concepts and the 
approximate theories could be checked. In the present 
paper, short of providing such an analytic solution, we 
obtain a very precise numerical solution to the problem 
of binary mixtures, by applying an extremely powerful 
cluster algorithm |8| to the problem. Like the algorithms 
of Swendsen and Wang || and Wolff |l(| for the case of 
the Ising model, the new method allows to obtain all the 
thermodynamic quantities with unprecedented accuracy. 
As a first important application, we actually establish 
the long sought for demixing transition both for spheres 
and for cubes. 

In our algorithm || , the convergence problems of ordi- 
nary Monte Carlo simulations are completely eliminated 
in a wide and physically interesting range of parameters. 
We have obtained convergence for up to 10 6 particles 
(limited by the size of the computer memory available 
to us) where previous work (for systems two orders of 
magnitude smaller) remained inconclusive. Our method 
applies equally well to spheres, cubes or any other shape, 
and both to the continuum and the lattice Jll[] . 

At any step of the algorithm (cf. j|] for details) one 
generates a new 'copy' of the 'original' by inverting the 
latter around a randomly chosen pivot point x p . Original 
and copy are then superimposed. The combined system 
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presents 'clusters' of overlapping objects. Some of these 
clusters are then 'flipped': particles belonging to the orig- 
inal are assigned to the copy, and vice versa. Thereafter, 
the copy is discarded and a new pivot point is chosen. 
The algorithm can easily be set up for simulations at 
constant particle number. It is completely symmetric 
with respect to the operation: original copy, a prop- 
erty which implies detailed balance. The method works 
perfectly well as long as the combined system breaks up 
into at least a few sizeable clusters. In our previous work 
on hard spheres in two dimension, we located this purely 
algorithmic percolation threshold at a much lower pack- 
ing fraction than would have been useful to study the 2- 
dimensional liquid-solid phase transition. In the present 
case of mixtures, the situation is vastly improved: we 
find that the percolation threshold (the optimal opera- 
tion point of our algorithm) mainly depends on the com- 
bined packing fraction r\ a + r]b, but not on the size ra- 
tio of the two species. This is radically different from 
the behavior of ordinary (local) Monte Carlo methods, 
where the diffusion of large particles becomes completely 
blocked by the nearby presence of many small ones. The 
lack of sensitivity of our algorithm to the size ratio of the 
particles is all the more interesting for a second reason: 
it has long been understood that the important phys- 
ical phenomena (depletion and the tendency to phase- 
separate) quickly move to lower overall packing fractions 
as the particles become dissimilar in size. These two ef- 
fects open up a large window of packing fractions and size 
ratios in which we can numerically solve the problem of 
liquid mixtures with the new algorithm. 

Let us first consider the superposed system of monodis- 
perse objects. In this case, the algorithm's optimal point 
of operation turns out to be at r\ a — r/p ~ 0.23 both for 
cubes and for spheres. The percolation threshold is thus 
located at a packing fraction corresponding to 1/4 of the 
close packing fraction for cubes and 1/3 for spheres. At 
?yp, the algorithm evenly flips clusters of any size. For 
larger packing fractions, intermediate cluster sizes will 
appear less often, since we either encounter the perco- 
lating cluster, or have to do with the algebraically de- 
creasing distribution of small clusters |12[ . Above the 
threshold, the algorithm deteriorates 'gracefully', and it 
is quite possible to converge the monodisperse system 
(at, say, N — 500) for packing fractions up to r\ ~ 0.4. 

We now introduce the very large number of small ob- 
jects (cubes or spheres). As long as the system remains 
homogeneous, we notice that the percolation threshold 
moves to slightly larger overall packing fractions. Dur- 
ing the simulation, we sample the partial distribution 
function of pairs of large particles gu(r) |13|. The nu- 
merical noise is much reduced if we consider not gu(r), 
but the integrated pair distribution function Gu(r) — 
inpi J r dr'r' 2 gu(r') (pi is the density of large particules). 
Gu (?") determines the average number of large particles 
within a radius r around a randomly chosen large parti- 



cle. We also compare gu(r) and Gu(r) to the pair dis- 
tribution functions g(r) and G(r) of the monodisperse 
system (with r/b = 0) at the same value of r\ a . Knowledge 
of the pair distribution functions for all r is equivalent to 
the computation of the structure factor (its Fourier trans- 
form), which informs us about the system's phase: for a 
homogeneous phase, gn(r) is completely flat for large ar- 
guments (Gu (r) will follow the monodisperse system's 
G(r)). In contrast, gu(r) and Gu(r) for phase separated 
systems will be system-size dependent functions even for 
moderate r. This fact translates to the presence of phase 
regions of varying extension. 
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FIG. 1. Integrated pair distribution function Gu(r) for 
spheres of radii r b /r a = 1/10 (N a = 62, N b = 62.000 and 
N a = 108, Nb = 108.000) (upper), compared to the monodis- 
perse case (lower). The packing fraction is rj a — nt = 0.1215, 
the system is homogeneous. 

For spheres at a packing fraction r\ a = rib = 0.1215, 
and a ratio of radii rb/r a = 1/10, we determined the in- 
tegrated pair distribution function Gu(r) and compared 
it to the monodisperse system (rjb = 0). In fig. 1, we 
present our results for N a — 62, Nb — 62.000 and for 
N a = 108, N b = 108.000 (the latter case, e. g., cor- 
responds to a box of side length L = 15.5 x r a , with 
periodic boundary conditions). We obtain very smooth 
curves, which indicate the exceptional convergence of the 
algorithm. More importantly, close agreement of Gu(r) 
with the monodisperse case for large r is reached. This 
clearly indicates that the introduction of small spheres 
has not changed the large-scale behavior of our system, 
which is still homogeneous. The inset of fig. 1 shows the 
Gu(r) for small arguments: the lower line corresponds 
to the monodisperse system, and the upper curve to the 
full simulation. A dramatic depletion effect is obvious. 
In agreement with previous knowledge, we observe that 



2 



the effect is strongest for values 2r a < r < 2(r a +r b ). Dif- 
ferentiating Gu (r) , we obtain gu (r) , whose contact value 
at r = 2r a is increased by a factor of 8.5 with respect to 
the monodisperse case, guir) oscillates for larger argu- 
ments, reaches a first minimum at r ~ 2{r a + r b ), and 
eventually levels out to the expected value gu(r) — » 1 
for large r. For small r, the integrated function, Gu(r), 
exceeds the monodisperse system's G(r) by about 0.4. 
In our opinion, this additional binding of an average 0.4 
particles characterizes the strength of the depletion much 
better than the contact value gu(2r a ). We also suggest 
that a system with an additional binding of more than 
one particle should be unstable to phase separation. 

The inset of fig. 1, as the main graph, contains in 
fact two sets of curves. On the scale of the figure, the 
curves for N a — 62 cannot be distinguished from those 
at N a = 108 (for r < L/2, neither can the data for 
N a = 864, which we have also computed). As mentioned 
before, the close agreement testifies to the good conver- 
gence of the algorithm, but also indicates that finite-size 
effects are completely negligible at these values of the 
physical parameters. This is a key observation, which 
leads us to strongly suspect that we are far away from 
any second-order phase transition point (as the critical 
point of phase separation) , which should lead to such ef- 
fects. 

Finally, at the combined packing fraction of rj = 0.243, 
we are still below the percolation threshold for the system 
with r b /r a = 1/10, but clearly above the monodisperse 
system's point of percolation. We believe that the lo- 
cal binding effects of the depletion force lead to a slight 
modification of the many-particle distribution functions, 
which is picked up in the distribution of cluster sizes. 

The situation appears to be changed for a ratio of 
radii r b /r a — 1/20. For this case, our memory re- 
sources have allowed us to perform calculations at N a = 
32, N b = 256.000 and N a = 62, N b = 496.000, again at 
Va = Vb = 0.1215 (L = 10.33 x r a and L = 12.88 x r a , 
respectively). Phase separation still does not seem to 
have taken place, but we notice the presence of important 
finite-size effects. For example, the system at N a = 32 
indicates an average 'additional binding' of 0.7 particles 
per sphere, while the larger system at N a = 62 yields a 
binding of 0.8. Again, the finite-size effects go into the di- 
rection of an increased depletion for the larger size. Since 
we are already at the limit of our computer resources, we 
were unable to check the phase behavior at even larger 
numbers of particles. 

Finally, we have performed simulations at r b /r a = 
1/30, and N a = 32, N b = 864.000, again at the same 
packing fraction rj a = i] b = 0.1215. In these very large 
simulations, the system clearly has crossed into the sepa- 
rated phase: the distribution function Gu (r) is suddenly 
very different from the monodisperse system's G(r) for 
all values of r, and the additional binding is clearly larger 
than 1. In this case, we usually observe the presence of 



one or two large aggregates |i~4j of large spheres which 
comprise most of them. In agreement with this obser- 
vation we also notice a dramatic change in the behav- 
ior of our algorithm: the distribution of 'cluster' sizes is 
shifted towards very large clusters, since the presence of 
the dense phase of spheres (the 'aggregate') pushes the 
system locally beyond the percolation threshold. 

Our algorithm is extremely powerful, but we have nev- 
ertheless reached the memory limits of today's worksta- 
tions. For spheres, the transition takes place at rather 
small size ratios (r b /r a ~ 1/20): therefore, we simulate a 
very large total number of particles but the Gu (r) belongs 
to a small system with only a few dozen large spheres. 
For spheres, the finite-size analysis, and the precise lo- 
cation of the critical point will have to be done on more 
powerful machines. 
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FIG. 2. G u {r) for cubes of side lengths d b /d a = 1/10 
(upper two curves at small r), compared to the monodis- 
perse case (lower) for N a = 62, N b = 62.000 and 
N a = 108, Nb — 108.000. The system, at a packing fraction 
of rj a = r) b = 0.067 is phase-separated. 

We have found it interesting to pursue our study in 
a different direction, i. e. to confirm phase separation 
in the model of hard parallel cubes of side lengths d a 
and d b , respectively. Here, the large cubes can touch 
on opposite faces. The osmotic pressure, exerted on a 
much larger surface, leads to a stronger depletion force, 
which should strongly favor phase separation. This is in- 
deed what we have found. As for hard spheres, we have 
observed i) depletion at any packing fraction; ii) finite- 
size effects, which become more important as the cubes 
grow more dissimilar in size, and which render the sys- 
tem more unstable for larger particle numbers; Hi) the 
phase-separated regime for very strong depletion. Again, 
the instability appears as soon as the average additional 
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'binding' is of one particle per cube. 

Simulations have revealed no signatures of an insta- 
bility for a size ratio of dt/d a — 1/2 (at the accessible 
packing fractions below r\ a + r/b = r/p ~ 0.23). Finite-size 
effects are negligible. For a ratio of db/d a — 1/10, at 
small packing fractions, the same holds true. However, 
for a packing fraction of r) a = rjb = 0.054, the system at 
N a = 62; L = 10.5 x d a remains clearly homogeneous, 
with a very strong depletion and an 'additional binding' 
of about 0.9 particles. At the same packing fraction, at 
N a = 496, N b = 496.000; L = 21.0 x d a , the G u (r) has 
pulled away from the monodisperse case for all r: the sys- 
tem has already undergone phase separation. In fig. 2, 
we present well-converged data for a slightly larger pack- 
ing fraction r\ a = rfb = 0.067, still at db/d a — 1/10. Here, 
already the system at N a = 62; L — 9.75 x d a is in the 
demixed phase. The data at N a = 108; L — 11.73 x d a 
illustrate the large finite-size effects at moderate r a . The 
maximum difference of (Ga — G) is much larger for 
N a = 108 than for N a — 62, and is expected to diverge for 
N a — * +00. We stress again that this effects are specific 
of the phase-separated system, and in the present form 
of the fixed particule number system, in which actual 
phase coexistence is obtained. In the demixed system, 
one of the two phases will be 'rich in cubes'. Unfortu- 
nately this phase will have a packing fraction above r/p 
(except very close to the critical point), and will be dif- 
ficult to study with our algorithm. The route towards 
instability is not touched by this problem. Nevertheless, 
both systems seem to have converged. While we did not 
study the demixed phases in detail, we have noticed the 
appearance of aggregates which are all but closed packed. 

In conclusion, we have studied the equilibrium phases 
of hard core mixtures. A superior algorithm has allowed 
us to establish and to access the demixed phase both 
for spheres and for cubes, and to investigate the subtle 
interplay between short-range depletion and long-range 
demixing. There are many questions and a large number 
of directions for further research, besides those already 
mentioned. Primarily, we think that the precise phase 
diagram needs to be established, especially the position 
of the critical point. In addition, the comparison with 
various closure formulas should be undertaken. We can 
already see that ref. Q places the critical packing frac- 
tion for cubes much too high. For most closure approx- 
imations, the numerical applications were done at quite 
large values of the size ratio, probably since Monte Carlo 
simulations for very dissimilar objects were thought to 
be completely out of reach. Paradoxically, the opposite 
is true. We are much more at ease at large asymme- 
try, as long as the packing fraction is not too high. The 
comparison between the exact numerical points and the 
closures should of course be done directly on the observ- 



ables, such as gu(r), and not on the phase diagram. To 
encourage and simplify further work, we will make avail- 
able via email the Fortran code used in this paper. 

Finally, the question remains whether the artificial per- 
colation threshold vp presents an unsurmountable barrier 
to the numerical solution. Several ideas to go much be- 
yond r]p have been formulated (cf Q). Even in the very 
dense limit of importance in the two-dimensional melt- 
ing problem, the flip of the percolating cluster can be 
avoided, but we have up to now been unable to trans- 
form this idea into a working algorithm Jl^] . However, 
we are firmly convinced that vp is not a hard boundary: 
many more difficult problems, such as melting, are also 
likely to fall under the attack of appropriate, while very 
specialized algorithms. 



buhot@physique.ens.fr; krauth@physique.ens.fr 
[1] S. Asakura and F. Oosawa J. Chem. Phys. 22 1255 
(1954). 

[2] Y. Mao, M. E. Cates and H. N. W. Lekkerkerker Physica 

A 222 10 (1995). 
[3] T. Biben, P. Bladon and D. Frenkel J. Phys. Condens 

Matter 8 10799 (1996). 
[4] J. L. Lebowitz and J. S. Rowlinson J. Chem. Phys. 41, 

133 (1964). 

[5] T. Biben and J. P. Hansen Phys. Rev. Lett. 66, 2215 
(1991). 

[6] M. Dijkstra, D. Frenkel and J. P. Hansen J. Chem. Phys. 

101, 3179 (1994); M. Dijkstra and D. Frenkel Phys. Rev. 

Lett. 72 298 (1994). 
[7] J. A. Cuesta Phys. Rev. Lett. 76 3742 (1996). 
[8] C. Dress and W. Krauth J. Phys. A: Math Gen 28, L597 

(1995). 

[9] R. H. Swendsen and J.-S. Wang Phys. Rev. Lett. 58 86 
(1987). 

[10] U. Wolff Phys. Rev. Lett. 62 361 (1989). 

[11] The first successful simulations using the algorithm were 
performed in a 3-dimensional lattice gas model: J. R. 
Heringa and H. W. J. Blote Physica 232A, 369 (1996); 
preprint to appear in Physica A; J. R. Heringa preprint. 

[12] D. Stauffer Phys. Rep. 54 1 (1979). 

[13] J. P. Hansen and I. R. Macdonald Theory of Simple Liq- 
uids, 2nd edition (Academic, London, 1986). 

[14] Two large spheres are 'aggregated' if their distance is 
smaller than 2(r a + r&). 

[15] S. Bagnier and W. Krauth, unpublished (1996). 

[16] The algorithm was carefully checked against a local 
Monte Carlo method for small systems. In one dimen- 
sion, excellent agreement with the exact solution |l7]] was 
obtained. 

[17] A. Robledo and J. S. Rowlinson Molec. Phys. 58, 711 
(1986). 



4 



